;;---------------------------------------------------------------------  
;; creat zonal mean radon and ipr distribution 
;;---------------------------------------------------------------------  
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" 
load "$KAILIBROOT/lib_interp.ncl"
load "$KAILIBROOT/lib_radon.ncl"


begin
  
  overland = True 

  path = "/pf/m/m222044/libdata/"
  pata = "/pf/m/m222044/rnpaper/"

  season = "JJA" 
  run = "RA007" 

  load "load_interp" 

  filenama = path+"mbase_T63L31.nc" 
  filenamb = pata+run+"_RN222_"+season+"avg.nc" 
  filenamc = path+"geop_T63L31_"+season+".nc"  
  filenamd = pata+run+"_RN222_"+season+"avg.nc" 
  filenamf = pata+run+"_RN222_"+season+"p90.nc" 
  filenamg = pata+run+"_RN222_"+season+"p50.nc" 

  fila = addfile(filenama,"r") 
  filb = addfile(filenamb,"r") 
  filc = addfile(filenamc,"r") 
  fild = addfile(filenamd,"r") 
  filf = addfile(filenamf,"r") 
  filg = addfile(filenamg,"r") 

  time     = fila->time
  lat      = fila->lat
  lon      = fila->lon
  phis     = fila->geosp
  slm      = fila->slm
  hyam     = filb->hyam
  hybm     = filb->hybm
  rns      = filb->RN222
  geop     = filc->var156
  ips      = fild->RN222
  rms      = filf->RN222 ;; 90pct 
  ros      = filg->RN222 ;; 50pct 

;;printVarSummary(slm) 
  print("--- data read ---") 


  rns = RN222_mBqm3(rns) 
  rms = RN222_mBqm3(rms) 
  ros = RN222_mBqm3(ros) 
  ips = RN222_mBqm3(ips) 

  print("--- conv unit ---") 


  ;; new height levels 

  hh = (/4000.,3500.,3000.,2750.,2500.,2250.,2000.,1750.,1500.,1250.,1000.,900.,800.,700.,600.,500.,400.,300.,250.,200.,150.,100.,60.,30.,0./) 
  
  nlat = dimsizes(lat) 
  nlon = dimsizes(lon) 
  ntim = dimsizes(time) 
  nlev = dimsizes(hh) 

  rnh = new((/ntim,nlev,nlat,nlon/),"float") 
  rnh = 0. 
  rnh@_FillValue = -999. 

  rnh!0 = "time" 
  rnh!1 = "lev" 
  rnh!2 = "lat" 
  rnh!3 = "lon" 

  rnh&time= time
  rnh&lev = hh
  rnh&lat = lat
  rnh&lon = lon

  rmh = rnh
  roh = rnh
  iph = rnh

  print("--- array created ---") 

  phisp = conform(geop,geop(:,30,:,:),(/0,2,3/)) 
 ;phisp = conform(geop,phis(:,:,:),   (/0,2,3/)) 

  ;; height above ground (m)  

  geop = geop - phisp + 60. ;;/ 10. 

  print("max,min of geop : "+max(geop)+"   "+min(geop)) 
  print("--- height caled ---") 

  ;;do m = 0, ntim - 1 

  ;;print("m = "+m) 

  ;; ion production rate 

  ips = RN222_ionprod(ips)  


  do j = 0, nlat-1
  do i = 0, nlon-1
      rnh(:,:,j,i) = (/ linint1( geop(:,:,j,i), rns(:,:,j,i), False, hh, 0 ) /) 
      rmh(:,:,j,i) = (/ linint1( geop(:,:,j,i), rms(:,:,j,i), False, hh, 0 ) /) 
      roh(:,:,j,i) = (/ linint1( geop(:,:,j,i), ros(:,:,j,i), False, hh, 0 ) /) 
      iph(:,:,j,i) = (/ linint1( geop(:,:,j,i), ips(:,:,j,i), False, hh, 0 ) /) 
  end do
  end do

  do k = 0, nlev-1  
     rnh(:,k,:,:) = where(ismissing(rnh(:,k,:,:)),rns(:,30,:,:),rnh(:,k,:,:))   
     rmh(:,k,:,:) = where(ismissing(rmh(:,k,:,:)),rms(:,30,:,:),rmh(:,k,:,:))   
     roh(:,k,:,:) = where(ismissing(roh(:,k,:,:)),ros(:,30,:,:),roh(:,k,:,:))   
     iph(:,k,:,:) = where(ismissing(iph(:,k,:,:)),ips(:,30,:,:),iph(:,k,:,:))   
  end do 

 
  ;;end do


  if (overland) then 

  do j = 0, nlat - 1
     if(all(slm(:,j,:).le.1.e-5))  then 
        slm(:,j,:) = 1. 
     end if  
  end do


  slmp = conform(rnh,slm,(/0,2,3/)) 

  ;printVarSummary(slmp) 
  ;printVarSummary(rnh) 

  rnh = where(slmp.lt.0.5,-999.,rnh) 
  rmh = where(slmp.lt.0.5,-999.,rmh) 
  roh = where(slmp.lt.0.5,-999.,roh) 
  iph = where(slmp.lt.0.5,-999.,iph) 

  end if ;; overland 


  print("--- output ---") 
   

 ;filenamo = "out_all_"+run+"_"+season+".nc" 
  filenamo = "out_"+run+"_"+season+".nc" 

  system("rm "+filenamo) 
  
  filo = addfile(filenamo,"c") 

  filo->RN222=rnh 
  filo->RN22290=rmh 
  filo->RN22250=roh 
  filo->IPRR=iph

end



